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ABSTRACT 

This paper describes a fourth-order integration algorithm for the gravita- 
tional A^-body problem based on discrete Lagrangian mechanics. When used 
with shared timesteps, the algorithm is momentum conserving and symplectic. 
We generalize the algorithm to handle individual time steps; this introduces fifth- 
order errors in angular momentum conservation and symplecticity. We show 
that using adaptive block power of two timesteps does not increase the error in 
symplecticity. In contrast to other high-order, symplectic, individual timestep, 
momentum-preserving algorithms, the algorithm takes only forward timesteps. 
We compare a code integrating an A^-body system using the algorithm with a 
direct-summation force calculation to standard stellar cluster simulation codes. 
We find that our algorithm has about 1.5 orders of magnitude better symplectic- 
ity and momentum conservation errors than standard algorithms for equivalent 
numbers of force evaluations and equivalent energy conservation errors. 

Subject headings: methods: N-body simulations — methods: numerical 



1. Introduction 

The gravitational iV-body problem is the numerical integration of trajectories for N 
particles with pairwise inverse square-law forces. Ideally, a numerical method should respect 
all of the symmetries of the exact problem. Conservation of linear and angular momentum 
and energy — in toto or in pairwise interactions — are often used as quality indicators for nu- 
merical algorithms. However, there is another conservation law following from Hamiltonian 
dynamics: conservation of the symplectic form on phase space. 



- 2 - 



Symplectic integration methods conserve exactly the symplectic form. This leads to 
other desirable properties of the flow of the integrator; for example, symplecticity implies 
incompressible, or dissipation less, flow. Symplectic methods a re well known in nunaerica l 
integration of the solar system following the pioneering work of IWisdom fc Holmanl (jl99ll ). 
However, the advantages of symplectic methods apply to systems with many more bodies. 
Indeed, cosmological simulations with millions or even billions of particles often use the 
leapf rog method because of its symplectic behavior (ISpringelll2005bl : IShirokov fc Bertschinger 
2005h . 



The leapfrog algorithm is second-order accurate. Higher-order accurate compositional 
algorithras are possible but co nventionally require sub-steps integrated backwards in time 
( lYoshida]ll993l : IChambersll2003l ) and such techniques typically require many more force eval- 
uations than standard integrators. Symplectic implicit Run ge-Kutta-Nystrom integrators 



are well-known, a nd involve relatively few force evaluations (jSurid Il989l : iMarsden fc West 
200ll : lstuchill2002h . 



Gravitating systems typically have a large dynamic range of density and hence dynamical 
time, making it computationally inefficient to use a constant timestep, as required by most 
symplectic algorithms. Either adaptive timesteps (which change with time as a system 
evolves), individual timesteps (which differ for each particle), or both are required to make a 
computation feasible. This is especially true when un softened inverse-squ are law forces are 
used, e.g., in numerical simulation of globular clusters (IHeggie &: Hutll2003l ). It is well-know n 
how to use compositional symplectic algorithms with individual timesteps (ISpringelll2005bl ). 
but the larger number of force evaluations in the high-order compositional algorithms and the 
necessity of back^yard sub -steps rule these out for use in simulations of large- iV gravitating 
systems (ISpringell l2005al ). No individual timestepping symplectic Runge-Kutta-Nystrom 
algorithm has appeared previously in the literature. 

Below we describe a fourth-order symplectic Runge-Kutta-Nystrom algorithm from a 
variational point of view. The algorithm requires the solution of a nonlinear algebraic equa- 
tion for one forwards sub-step. If the nonlinear equation is solved approximately by iteration 
then the algorithm is approximately symplectic with an error in phase space conservation 
that can be made arbitrarily small. We describe how to use data from the last step of the 
integrator to generate an approximate solution sufficient to obtain fifth-order symplecticity 
and momentum-conservation using only two force evaluations per step. We generalize this al- 
gorithm to use individual (and adaptive) timesteps; the formulation in terms of a discretized 
action principle is essential to the generalization. This is the first individual timestepping 
symplectic Runge-Kutta-Nystrom algorithm to appear in the astrophysical literature. 



It is widely believed that adaptive timesteps are incompatible with symplecticity but 
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we show otherwise in Section HI This assumed breakdown of symplecticity, coupled with 
the large number of force evaluations for standard higher-order compositional symplectic 
methods and the lack of individual timestepping symplectic Runge-Kutta-Nystrom algo- 
rithms, has led resear chers to seek alternative ways to reduce diss ipation, e.g. using time- 



reversible integration (IMakino et al.l Il996l : iPreto fc Tremaind Il999l ) . These techniques treat 



the symptoms of non-symplecticity (linear growth in various errors) without treating the 
cause (non-conservation of the symplectic form). Our algorithm addresses the cause, and we 
see corresponding improvements in symplecticity and momentum conservation for equivalent 
energy error to relative standard algorithms. 

In this paper we present a fourth-order integrator requiring only two force evaluations 
per timestep, which is fifth-order in symplecticity. Each additional force evaluation improves 
the symplecticity by two powers of the timestep. We generalize the integrator to individ- 
ual timesteps and analyze the breakdown of symplecticity when adaptive and individual 
timesteps are used. We show that symplecticity is effectively restored when block power of 
two timesteps are used. 

The algorithm is based on a discrete approximation to the action of a system, described 
in Section [2l The algorithm contains a non-linear equation whose solution must be approx- 
imated; we compare two approximation methods in Section [31 Adaptive, individual, and 
combined block timesteps are discussed in Sections HI O and [6] respectively. Numerical tests 
are presented in Section [71 Conclusions are given in Section [HI 



2. Variational Integrators 



Variational integrators are based on applying Hamilton' s principle o 



discrete approximations to the action for a physical system. [Lew et al.l (|2004|) is an excellent 



stationary action to 



introduction to variational integrators in an engineering context; iMarsden fc WestI (l200ll ) 
provides a much more mathematical discussion, including proofs of the essential properties 
of variational integrators and many examples of particular integration rules. This section is a 
brief introduction to variational integrators. Here and throughout we suppress vector indices 
on variables (juxtaposition of variables thus denotes multiplication in the one- dimensional 
case and the usual dot-product in the multidimensional case). We denote the derivative of 
the function / by Df; we denote the partial derivative on the ith argument of the function 
g by diQ (argument labels begin at 0). 



The fundamental theorem of variational integration (IMarsden &: West 



2001 



, Theorem 
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2.3.1) states that if H is an approximation to the action of a mechanical system, 
H{h, gi, g', g", . . . , g("), q^) = S[q]{UM + h)+0{K'-+^) = T^' rft g(t), Dg(t)) + 

(1) 

where gi = g(ti), g2 = g(^2)5 the g^*) are intermediate positions in the time interval + /;,], 
S is the action functional and L is the Lagrangian for the system, then the equations 

9iif(/i,gi,g',g",...,g("),g2) = -pi (2a) 
d,H{h,q^,q\q'\...,q^^\q2) = 0, z = 2, 3, . . . , n + 1 (2b) 

a„+2i^(/i,gi,g',g",...,g("\g2) = P2 (2c) 

define a map (gi,pi) ^-^ (^'2,^2) which is an order r integrator for the mechanical system. 
The function H is called the discrete Lagrangian. Equations (12bp extremize the discrete 
action approximation with respect to the discrete path {g*-*^}, while equations fl2aj) and fl2cl) 
exploit that the action is a Fi-type generating function for the time-evolution canonical 
transformation (jSussman et al.ll200ll . pp. 415-416). 

The map defined by equations ([2]) has many useful properties analogous to the properties 
of the exact evolution of the system defined by L. First, it is momentum-preserving: imagine 
an infinitesimal variation in the coordinates gi, g2 and g'-*-' which leaves H invariant when 
gi, the g*^*^ and g2 satisfy equations ([2]). We have 

n+l 

5H = d,H{h, gi, g', g", . . . , g("), g2)5gi + ^ diH{h, gi, g', g", . . . , g^'^), g2)5g« 

i=2 

+ dn+2H{h, gi, g', g", . . . , g^"), g2)5g2 
= P2Sq2 - PiSqi 

= 0. (3) 

In this situation the quantity p5q is conserved by the integrator: this is the discrete version 
of Nother's theorem. Assuming that H inherits the symmetries of L, the integrator will 
exactly preserve the associated discrete momenta. 

Second, the map is symplectic. Consider the discrete approximation to the action over 
an interval evaluated on the integrator path: 

M-l 

5 (gi, g2, . . . , gjw) = (h, qi, q-, . . . , gf \ g^+i j , (4) 

i=l 

where the g^ satisfy the integrator equations ([2]). Taking one exterior derivative of S gives 



c^5(gi,g2,...,gM) = d,H (h, q,, q[, . . . , q^^\ q^) dq. 
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(5) 



the terms involving dq2, . . . , dqu-i are identically zero because the trajectory satisfies equa- 
tions fl2aj) . pbj) and fl2c|l . Two exterior derivatives of S give zero, yielding 

d'S{quq2,....qM) = dA+^H [h, q,, q[, . . . , q^^\ q^) dq, A dq^ 

+ didn+2H (ji, qu-i, q'u-i^ • • • > 9m-1' 9m j <igM A dq^ 



M-l 



0. 



(6) 



Instead of considering evolution on phase space consider the corresponding evo- 

lution on the discrete state space {qi,q2)- Evolution maps the initial state-space for H, 
(gi, ^2) £ IR™ X IR™", to an isomorphic space, (^m-i, Q'm) G x M™', where m is the dimen- 
sionality of configuration space. Equation can be written 

didn+2H (h, qi, q[, . . . , qf\ 52) dqi A dq2 = F* didn+2H (h,qi,q[,..., q^\ ) dqi A dq2 



(7) 

using the pushforward map under evolution, F* . All forms in equation ([7]) live on the 
cotangent bundle of the state space M*" x R™. We see that the integrator conserves the 
discrete symplectic form on the state space of H, 



dA+2H (/.,gi,g;,...,gS"\g2) dq^ A dq2. 



This is the direct analog of the symplecticity of continuous time-evolution in a Hamiltonian 
system. 



Using equation fl2a|) . we see that 



- dpi = didiH{h, gi, q', q", . . . , q^^\ q2)dqi + dn+2diH{h, gi, q', q" , . . . , g2)c^g2, (9) 

and therefore conservation of the discrete symplectic form in equation ([8]) implies conserva- 
tion of the Poincare integral invariant on phase space: 



(n) 



didn+2H {h, qi, q[, q^^\ ^2) dqi A dq2 = dpi A dqi. 



(10) 



Finally, while it is impossible for a constant-time-ste pping integrator to be momentum- 
conserving, symplectic and to exactly conserve en ergv (IGe fc MarsdenI Il988l ). variational 
integrators generally have bounded energy error. iLew et al.l (120041 ) explain that the dis- 
cretized trajectory is sampling the continuous trajectory of a Lagrangian system, L, which 
is near L. L satisfies 

dt L{t, q{t), Dq{t)) = H {h, qi, q' , q" , . . . , ^2) , (H) 
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where the integral is evaluated on the trajectory which satisfies the Euler-Lagrange equations 
for L with q{ti) = qi and q [ti + h) = q2 and the positions which are arguments for H 
satisfy equations ([2]). Equation ([TT]) implies that H is the exact generating function for time 
evolution under L. In general it is only possible to compute a truncation of L to any desired 
order in h, but it is possible to prove that L is close to L in the space of possible Lagrangians. 
Since the trajectory remains on the energy level set of L in phase space, near the energy 
level-set for L, the energy error remains bounded. 



2.1. Galerkin Gauss-Lobatto Variational Integrators 



To define a variational integrator, we need an H-type function which approximates the 
action for a syst e m ove r an interval. There are many ways to find such a function; see 
Marsden fc WestI (120011 ) for an extensive discussion of the various types of integrator. In 
this paper, we will focus on the so-called Galerkin Gauss-Lobatto (GGL) integrators. These 
integrators assume a polynomial trajectory in time and approximate the action integral using 
a Gauss-Lobatto quadrature rule. Gauss-Lobatto quadrature is appropriate because it gives 
the highest-order integration rule for a given number of points subject to the constraint 
that the Lagrangian is evaluated once at the beginning and once at the end of the interval. 
Evaluating at the beginning and end of the interval is important because it preserves the 
symmetries of the continuous Lagrangian. 



Marsden fc WestI (120011 ) show that all GGL integrators can be written as symplectic 



partitioned Runge-Kutta integrators and derive formulas which relate the SPRK coefficients 
to the discrete Lagrangian. Because variational integrators are symplectic, GGL integra- 
ly sa ti sfy the constraint s on s y mplecti c Run ge-Kutta coefficients (see, for 



tors automatical 



example, ISurij (119891 ): iMarsden fc WestI (120011 ): IStuchil (120021 )). The particular integration 



algorithms we consider here belong to a sub-class of partitioned Runge-Kutta methods often 
called Runge-Kutta-Nystrom methods. 

We define 

Hih,q,,q',...,q^-\q2)^ 

wiLiti, gi, g', . . . , q^'^\ q2, h.h), do(f){ti, gi, g', . . . , g^"\ g2, ti, h)) 

n 

+ w^U{t^^, gi, g', . . . , g("), g^, h, h), 9o0(t», gi, g', . . . , g("), g^, h, h)) 

i=l 

+ W2L{ti + h, (j){ti + h, gi, g', . . . , g^"), ga, h, h), do(j){ti + /i, gi, g', . . . , g("), ga, h, /i)](,12) 

where </)(t, gi, g', . . . , g^"-*, g2, h, h) is the interpolating polynomial for q{t) passing through the 
points {gi, g', . . . , g("\ g2} at times {ti, t', . . . , t^'^\ti + h}. The times {ti, t', . . . , t^'^\ti + h} 
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and weights {wi, w', 



w 



in) 



W2} define the n+2 point Gauss-Lobatto quadrature rule. This 



rule is exact for quadrature of polynomials up to and in cluding order 2n + 2. The inte grator 
so defined will be of order 2n + 2. See, for example, lAbramowitz fc StegunI (Il972l . Table 
25.6) for appropriate times and weights. This is easier than it looks; examples follow. 



2.1.1. Two-Point Integrator 



The two-point Gauss-Lobatto integration rule has times {ti, ti+Zi} and weights {h/2, h/2}] 
a two-point interpolating polynomial is a line. Therefore, we have 

h ' 



H{h,qi,q2) = 

For a Lagrangian L(t,q,v) 
integration rule 



L ti,qi 



q2 - qi 



+ L [ ti + h, q2. 



q2 - qi 



(13) 



h J ■ " ■ h 
|mt>^ — V{q), equations (l2al) and (l2cl) result in the explicit 



?2 



P2 



qi + h 



Pi 



m 



2m 



DV{q, 



h 



p,--[DV{qr) + DV{q2)]. 



(14a) 
(14b) 



This rule is kick-drift-kick leapfrog; it can be derived from the Hamiltonian viewpoint by 
iterating the evolutions of the splitting of the Hamiltonian for this system into Hi{q,p) = 
V(q) and H^jq^p) = p^/ (2m) (for a thoroug h exploration of this idea see, for example. 
Wisdom &: HolmarJ (jl99ll ) and lYoshidal (119931 )). This rule is second order, as expected from 
the order of the quadrature rule. 



2.1.2. Three-Point Integrator 



The three-point Gauss-Lobatto quadrature rule has times + h/2,ti + h} and 

weights {h/Q,2h/3, h/Q}. The three-point interpolation polynomial is quadratic in time. 
Applying equation ( |T2l) . we find 



H{h,qi,q',q2) = h 



L ti,qi 



Aq' - 3qi - q2 



h 



+ -L ti + -,q' 



h 

g2 - gi 
h 



-L{ti + h,q2 



qi + 3^2 - 4g' 



h 



liL{t,q,v) 



V{q) then equations ( l2al) . ( l2bl) and ( l2cl) reduce to 



qi + 



hpi 1 f h 
2 m ~ I 



2DV{qi) , DV{q') 



3m 



3m 



.(15) 



(16a) 
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92 

P2 



m I 



DV{qi) 2DV{q') 



3m 



3m 



Pi 



h 



DVjqi) ^ 2DV{q') ^ DV jq,) 



6 



6 



(16b) 

(16c) 



To implement this integration scheme, we must solve the implicit equation fll6al) for q'. In 
spaces of high dimensionality (i.e. an A^-body system with large A^), the only efficient way 
to do this is through iteration: we treat equation fll6ap as a prescription for updating the 
value of q', and iterate. We describe two iteration techniques in the following section. 

The three-point GGL integrators take only forward steps in time. It is not possible to 
write these integrators as iterated evolution of a splitting of the Hamiltonian corresponding to 
L (as in the two-point case): a general theorem states that it is impossible to have evolution 
by Ha miltonian sp l itting which is higher than second-order accurate and takes only forward 
steps JsheneJll989l : ISuzukil EoOll ) . 



It is possible to formulate higher-order mapping integrators fr om the Hamiltonian 



perspective which take only forward st eps using the force gradient (IWisdom et al.l Il996 



Scuro fc Chin 



20051 : IChin fc Chenll2005l ). but force gradients can be expensive to compute. 
Omelyan J2006h describes how to approximate computation of a force gradient with an ex- 



tra force evalu ation at a s hifted position. The iteration technique in Section [3?T] reproduces 
algorithm 8 in lOmelyanI (120061 ) from the equations fll6ap . fll6bl) . and fll6cp . However, we 
find that for equivalent numbers of force evaluations in the A^-body problem, the alternate 
iteration technique in Section 13.21 typically outperforms the one in Section 13.11 by one to two 
orde rs of magn i tude in energy error (see Figure [1]) , so the gradient-approximating algorithm 
from lOmelyaru (120061 ) is sub-optimal for our purposes. 



Forward time steps are important for cosmological simulations which include gas dy- 
namics because such simulations are unstable under time reversal. The requirement of 
symplecticity and forward timeste ps has previous ly restricted cosmological simulations to 
second-order mapping integrators (jSpringell l2005al ) . 



It is straightforward to derive the integration equations for n + 2-point GGL integrators 
with n > 1. All such integrators have implicit equations for the intermediate positions, 
g', . . . , q^'^\ which must be solved via iteration exactly as in the n = 1 case discussed above. In 
the presence of individual timesteps (Section[5]), we must predict the intermediate positions — 
we cannot iterate the implicit equations to convergence. For integrators of order greater 
than four, such predicted positions do not solve the implicit equations accurately enough to 
make the symplecticity error scale better than the trajectory error; as we shall see, we can 
predict the solution to equation fll6ap accurately enough to produce fifth-order symplecticity 
error in the fourth-order integrator. For this reason, the fourth-order integrator is uniquely 
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positioned in the hierarchy of GGL integrators for adaptation to individual timesteps 



3. Solving the Implicit Equation 



This section discusses two ways to solve equation fll6al) with iteration. They differ in 
their choice of initial guess for q' . The choice in Secti on [3.11 produces an algorithm which 
is compositional, and is equivalent to algorithm 8 from lOmelyanI (120061 ). That algorithm is 
exactly phase-space-volume and momentum conserving, but only fourth-order symplectic, 
with three force evaluations per step. The choice in Section [3l2] produces an algorithm which 
is fifth-order symplectic (and the same in phase-space- volume error), exactly conserves linear 
momentum, and conserves angular momentum at fifth-order, with two force evaluations per 
step. 



OmelyanI (120061 ) reports excellent energy conservation for the algorithm in Section 13.11 



when simulating the one-dimensional Kepler problem (where phase-space-volume conserva- 
tion implies symplecticity) , but we find that the algorithm in Section [3.21 has one to two 
orders of magnitude better energy conservation for equivalent numbers of force evaluations 
in the iV-body problem for N > 2 (see Figure [1]). This is probably due to the superior 
symplecticity of the algorithm in Section 13.21 for multidimensional configuration spaces. We 
do not discuss the generalization of the algorithm in Section 13.11 to individual and adaptive 
timesteps, but instead focus on the algorithm in Section 13.21 for the remainder of the paper. 



3.1. Compositional Algorithm 



To reproduce algorithm 8 in lOmelyanI (120061 ). let the initial guess for q' be 



hpi 1 fhy2DV{qi] 



2m 2 

Then iterate equation (I16al) once, producing 



3m 



(17) 



Q = (1(0) 



24m 



DV (g[o)) • 



(1^ 



q2 and p2 are then given by equations (]16bl) and (I16cl) . with this choice for q'. This choice 
of initial guess and single iteration allows the algorithm to be written compositionally as 



P 



p-^DV{q) 
6 



(19a) 
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h p 

2 m 
2h 



P 



P 



DV (q- 



24m 



DV{q) 



P 



h p 

/ m 

p-^DV{q). 
b 



(19b) 
(19c) 

(19d) 
(19e) 



This is exactly the sequence of operations inlOmelvanI (120061 ). algorithm 8; equation (I19cl) 
is the approximation derived in Omelyan ( 20061 ) to the force gradient required in the corre- 
sponding algorithm of I Chin fc ChenI (120051 ). The algorithm needs four force evaluations for 
a single step, but only three in a long-running simulation because the first force evaluation 
of a step occurs at the same position as the last force evaluation of the previous step. The 
algorithm exactly conserves phase-space volume and momentum, but is only fourth-order 
symplectic (in 2n-dimensional phase-spaces with n > 1). As stated above, in practice we 
find that the energy error from this algorithm in A^-body simulations is significantly worse 
than that from the following algorithm. We compare the energy error behavior of the algo- 
rithms in Figure [H 



3.2. Prediction Algorithm 

As we shall see in Section |5l it is not, in general, possible to iterate equation (I16al) in the 
presence of individual timesteps. Iterating equation (I16ap for a particle corresponds to re- 
running the evolution of that particle over the first half of its timestep. This may require re- 
running many entire steps of particles with smaller timesteps than the given particle, which, 
in turn, may require re-running yet more steps of particles with even smaller timesteps. The 
explosion of work is exponential in the number of distinct timesteps assigned to particles in 
the simulation. 

Given that we cannot iterate equation (I16al) . it makes sense to try to use all the infor- 
mation at hand to predict the solution q' as well as possible at the beginning of a particle's 
timestep. The initial guess 

where t is the time at which the particle is at position qi, and F{t) is the force on the particle 
as a function of time, solves equation (I16al) with an error term of order h^. (Note that the 
final term in equation ( !20l) is twice what would be expected in a power series for q{t + h/2).) 
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We can use the force evaluations from the previous step, at times t — h, t — h/2 and t, to 
approximate F{t), DF{t) and D'^F{t) to sufficient accuracy to compute equation fl20|) with 
no new force evaluations. 



The integration algorithm 



92 



P2 



qi + h— 
m 



2 



Pi-h 



DVjqi 
6 



DVjqi 
3m 



2DV 



2DV 



+ 



3m 



DViq2) 



6 



(21) 



(22) 



which is equations ( |T6l) with q' replaced by the predicted g^Q-j from equation ( l20l) . is fifth-order 
accurate in symplecticity and angular momentum conservation, exactly conserves linear mo- 
mentum, and requires only two potential evaluations per step in a long-running simulation. 
It outperforms the algorithm in Section 13.11 in energy error for equivalent force evaluations 
by one to two orders of magnitude, as evidenced in Figure [H This is the algorithm we shall 
discuss for the remainder of the paper. 



4. Adaptive Timesteps 

In A^-body simulations it is essential to be able to adapt the timestep taken by an 
integrator to the local conditions of the system. Formally, we can model such an adaptive 
step by treating the parameter h in the discrete Lagrangian describing an integrator as a 
function of the positions through the step: 

H (h, gi, q[,..., q[^\ q^^ ^ H {h [qi, q[, . . . , q["'\ q-^ , gi, gj, . . . , gi"-", g2) • (23) 

Because the sequence of positions gi, gj, . . . , q^\ q2 encodes all available information about 
the trajectory, essentially any technique for choosing a timestep can be recast in this fashion. 
We do not change the integrator equations 

aii/(/z(gi,g;,...,g('^),g2),gi,gl,...,gi"\g2) = -pi (24) 

9,/7(/i(gi,g;,...,g("),g2),gi,g;,...,gi"\g2) = 0, z = 2, 3, . . . , n + 1 (25) 

dn+2H (h{qi,q[,...,q^''\q2) ,qi,q[,...,q^^\q2^ = P2- (26) 

Assuming that the timestep function h ^gi, q[, . . . , q^\ g2^ is invariant under the same 
coordinate variations which leave H invariant (this will be the case if both H and h inherit the 



Energy error vs force evaluations for different iteration tecliniques 




10 10 
Force Evaluations 



10 



Fig. 1. — Energy error versus number of force evaluations for the algorithms described 
in Sections 13.11 and 13.21 in a simulation of a 100-body Plummer model initial condition 
with constant, shared timesteps over a total time interval T = 1.0. As described in the 
text, the algorithm in Section 13.21 outperforms the compositional algorithm in Section 13.11 
for equivalent numbers of potential evaluations, though both algorithms exhibit the same 
asymptotic scaling. 
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same symmetries from the continuous Lagrangian), then the proof of momentum conserva- 
tion in equation ([3]) still holds. All terms in 5H proportional to OqII {h ^gi, q[, . . . , , Q'l, 
vanish because 6h also vanishes. Adaptive stepping poses no threat to momentum conser- 
vation, provided that the steps are chosen in a way which respects the symmetries in the 
continuous mechanical problem. 

Unfortunately, adaptive stepping does pose a threat to symplecticity. With a timestep 
function h (^qi, q[, . . . , q^\ ^2^, the expression for dS (see eq. ([5])) gets additional terms: 



M-l 



1=1 
X 



J2doH{h{q.,qi,...,q^\q.^?l,q.,q',...,qT\q. 

^doh (^qi, q'i, g^+i j dqi + YTj=i 9jh (qi, q-, . . . , gi+i) rfgp^ 

+ dn+ih (^qi, g-, . . . , g^.+i) dq,+i . (27) 

These extra terms prevent us from writing d'^S = as the difference of two forms, one of 
which is the pushforward of the other, as we did in equation (jTj). With general adaptive 
timesteps, there is no two-form on the state space which is conserved over a fixed number of 
steps. 

This can be understood intuitively in the following way. The conservation of the 
Poincare integral invariant, 

I = ^dq' Adpi, (28) 

i 

with the sum taken over coordinate dimensions, along any trajectory implies symplecticity 
(recall that eq. ffTOj) shows that conservation of the two-form didn+2H {h, qi,q[, . . . , qf \ g2 j dqi/\ 
dq2 on state space implies conservation of I on phase space). Note that rfg* and dpi are ex- 
terior derivatives and that X is a two-form. The Poincare integral invariant measures the 
sum of the areas of a tube of trajectories infinitesimally near a reference trajectory projected 
onto the sub-phase-spaces {q\pi)- But an integrator with an adaptive timestep does not 
advance all trajectories in the tube with the same h. Even if we had an adaptive timestep 
integrator which implemented the exact, continuous evolution of the system it still would 
not conserve I over a single step — stopping the evolution of the different trajectories in the 
tube at different times spoils the symplecticity of the continuous evolution. For an adap- 
tive timestep integrator symplecticity after a fixed number of steps is the wrong condition 
to consider. Rather, we should ask whether the integrator conserves the symplectic form 
over a fixed total time of evolution. In general, any well-behaveci^ integrator (including a 



^In this case, "well-behaved" implies an integrator whose maps converge uniformly at order to the 
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variational integrator with general adaptive stepsizes) which has trajectory error of order h"^ 
conserves the symplectic form at least to order h"^ in this sense. 

For variational integrators which choose timesteps using the popular block-power-of- 
two scheme we can do better: these integrators conserve the symplectic form almost every- 
where in phase space. In the block-power-of-two scheme, a function h^nax (^Qi, Qi, ■ ■ ■ , <ii'\ Q-i^ 
limits the maximum timestep; the actual timestep taken is rounded down from /imax to 
the nearest number of the form T/2", with n an integer, and T some total evolution 
interval. If the function /imax is continuous, about every point in state space for which 
^max [qi^ Qi, ■ ■ ■ , 7^ T/2" for some n there is an open neighborhood of points which 

round down to the same timestep. Thus, the actual timestep function h (^qi, q[, . . . , 
is piecewise-constant on state space, and the derivatives in equation fl7r|) vanish almost 
everywhere on state space for each step. A variational integrator with adaptive block-power- 
of-two timesteps is therefore symplectic almost everywhere on state space because it is a 
composition of symplectic steps. ("Almost everywhere" should be taken in the mathemati- 
cal sense of "except on a set of measure zero.") In Section [7] we present numerical evidence of 
this phenomenon. Incidentally, block-power-of-two timesteps have advantages for paralleliza- 
tion which have led other auth ors to adapt time-symmetric methods to block-power-of-two 



timesteps (iMakino et al.lll996l ). Variational integrators fit naturally into a block-power-of- 
two scheme which, in addition to its ease of parallelization, preserves symplecticity in the 
presence of adaptive timesteps. 



5. Individual Time Steps 

Astrophysical simulations of realistic systems often require individual time steps for each 
body or for subgroups of bodies. A variation of the three-point GGL integrator described 
in Section 13.21 allows for individual steps. This appears to be the first individual timestep 
Runge-Kutta-Nystrom integrator to appear in the literature. 

Our generalization to individual timesteps will exploit the way in which we derived the 
3-point GGL integrator. In the derivation, we assumed a quadratic-in-time trajectory for 
the particle which interpolates between the positions qi at ti, g' at ti + /i/2, and q2 at ti + h. 
We will write 

g(t) =0(t;gi,g',g2,ti,/i) (29) 
to represent this interpolated trajectory in what follows. We now repeat the derivation which 



exact evolution map and whose derivatives converge uniformly as well. 



- 15 - 



led to equations f[T6|) . allowing each body to have a different timestep. 

Write the iV-body Lagrangian as a sum over particle indices i = 1 . . . N: 



N 

E 

i=l 



N 



j=i+i 



m'm^V {q\t)-q^{t)) 



N 



(30) 



j=i 



where Visi) = "l/kl (recall we suppress vector indices). Assume for the moment that the 
timesteps of each body, h\ are fixed in time (Section [6] discusses how to integrate adaptive 
steps into this algorithm), and that particles are sorted so that < when i < j. Note 
that T* depends only on the trajectory of body i, while V* depends on the trajectories of 
bodies j with j > i. 

Now compute the discrete approximation to the action over a total time interval At 
using the three-point Gauss-Lobatto quadrature rule to integrate each of the — V*: 



N n'-l 



i=l k=0 



6 



+ 



2 {T' - V' 



t=kh' 



+ 



6 



(31) 



t={k+l)hi 



where At = n'^h'^ for each /i*. Recall that V* contains contributions from potentials V {q'^{t) — q^{t)) 
with i > i] this sequence of potential evaluations is illustrated in Figured 

Discretize the trajectories of each particle using a quadratic interpolation between the 
integration points: 



q\t) = {t; ql qt ql^„ kh\ h') , kh' < t < {k + l)h\ 



(32) 



A given interpolation point, ql, g^* or ql^i, will appear in T* and V* on step k of index i, 
and V\ j < i, on all steps / such that Ih^ < (A; + ^)h^ and {I + l)h^ > kh^. In general, ql, q'j^ 
or ql_^_-j^ appear together in as part of an interpolation 0, while each appears once alone 
in the three V* terms of timestep k. Figure [3] illustrates graphically the time sequence of all 
potential evaluations involving three particular interpolation points, ql, c(l and g^+i- 

Extremizing the action with respect to g^, g^* and g^^_]^ as in Section [2], and using equation 
to predict g^* as 



4 



gives 



i h'pl 1 



'''^\{kh^) + U^yDF{kh^) 



1 

12 



D^F {kh') , (33) 



<tk + ^ 



Pk 

m 



{hi' 



3m I dql 



E 

j<i 



dql 



3m \ dq 



j<i 



(34a) 
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Time 



Particle Index 



Fig. 2. — Time-sequence of potential evaluations in the integral approximation to the action, 
equation (!3T|) . assuming, for graphical simplicity, that = 2^~^h^. Time runs vertically and 
particle index runs horizontally. A right-facing arrow crossing a particle timeline represents 
evaluation of potentials involving that particle and all particles of greater index. Long 
arrows represent evaluations of potentials at the beginning and end of particle timesteps 
(associated with coefficients h/6)] short arrows represent evaluations of potentials at particle 
half-timesteps (with coefficient 2h/3). 



Time 



Particle Index 

Fig. 3. — Graphical representation of the potential evaluations involving three particular 
interpolation points, ql, q'^ and for particle i. Time flows vertically, and particle index 
runs horizontally. Potential evaluations in V go to the right; each of these evaluations 
involves exactly one of the positions g^, q'j^ and q^j^i- Potential evaluations in , j < i, come 
from the left, and generally involve the interpolated position of particle i. 
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Pk+l 



Pk 



+ 



6m 
dV 



3m \ dq'i^ ^ dq'^ 



j<i 



6m \ dq 



fc+i 



j<i 



(34b) 



We can evaluate the force terms in the above equations in the same time sequence as the 
potential evaluations in equation (|3T|) . First, assume we have predicted ql^^ using F^kh^), 
DF {kW) and D^F (kh^) in addition to predicting g^*. As we advance the system forward in 
time (imagine a horizontal line rising upward in Figure[2]), we encounter potential evaluations 
between various pairs of particles i and j. When we encounter a potential evaluation, we 
must compute the corresponding contribution to the dV/dq terms in equations We 
use the chain rule: each dV/dq term is the product of a force evaluated at interpolated 
positions and the interpolation coefficient of the corresponding coordinate in the position 
interpolation. We accumulate all such contributions over the course of each timestep. 



When we reach the end of a particular particle's step, we use equations fl34j) to advance 
the position and momentum of that particle, and then predict the new 1k+2 using 

the stored force data. Each particle must store its mass, its timestep, h\ three positions 
for the timestep, ql, q'^ and g^+i, the initial momentum and three force accumulators, 
dV/dql, dV/dq'^ and 9V/(9g^^^, for a total storage of 23 floating-point numbers in three 
dimensions (in double precision this is 184 bytes per particle). 



6. Adaptive Stepping and Individual Block Power of Two Timesteps 

In this section we combine the block-power-of-two adaptive timesteps described in Sec- 
tion H] with the individual timestep algorithm described in Section O In order that equation 
( |3T|) approximate the action for the A^-body system, we require /i* < for i < j. Thus, 
changing the timestep for a body potentially requires re-ordering the sum in equation (130|) . 
Bodies i and j can only change relative position in the sum when both have flnished a step, or 
there will be some potentials whose evaluations in the discrete Lagrangian do not correspond 
to a proper Gauss-Lobatto quadrature rule. Using block-power-of-two timesteps. 

At 

h' = i^. p' = lor2or3..., (35) 

for some total integration interval At. With block-power-of-two timesteps, all bodies com- 
pleting a timestep at a given time are contiguous in the sum from indices 1 to some We 
can re-compute the maximum timestep allowable for each of these bodies and sort them in 
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the sum according to their maximum timestep. Once the maximum timesteps are calculated, 
we must ensure that the actual timesteps taken satisfy 

h'<mm{hl^,^,h^'^'''-''), z = l,2,...,z-- (36) 



subject to the power-of-two restrictions in equation fl35|) . This procedure preserves the 
invariant that < for i < j and ensures < h\^^^. 



Makind (jl99ll ) recommends choosing a timestep which reflects the error of the predictor 



relative to the final solution. We choose 



/^Lx = I ?r'f ' (37) 

where t] is an accuracy parameter of the simulation, and gpj.^^ is the predicted position of 
the particle using equation fl2U]) . while is the position as determined by the integrator at 
the end of the prior timestep. The exponent is 1/5 because the predictor has an 0{h^) error 
term. 

The complete algorithm to advance body i by dt, assuming an array of bodies sorted 
by /imax and indexed from to — 1 is then: 



1. l{dt< hi^^^ then 

(a) Predict the positions of body i at t + dt/2 and t + dt using equation (120!) . and a 
Taylor series for q(t + dt), respectively. 

(b) Compute potential gradients between body i and bodies j, with j > i, at times 
t, t + dt/2, and t + dt, distributing the forces across the dV/dq accumulators of 
body j according to the interpolation coefficients for the position of body j at 
these times. Also accumulate the forces into the corresponding accumulator of 
body i. 

(c) Unless i = 0, advance body i — 1 by dt. Upon completion of this step the forces 
on body i from all bodies with indices below i will be stored in i's accumulators. 

(d) Update the position of body i using equations (IMI) . Calculate a new /i^ax- 

2. Otherwise 

(a) Advance body i by dt/2. 

(b) Sort bodies to i according to their maximum timesteps. 



(c) Advance the body at index i (which may not be the same body after the sort) by 
dt/2. 



The algorithm begins by attempting to advance body — 1 by the entire time interval, At. 
The order of computations and recursion above ensures that the timesteps are block-power- 
of-two, that the total integration terminates after advancing the bodies by exactly At, and 
maintains the invariant that < when i < j. 



In this section we report on numerical experiments involving two different GGL varia- 
tional integrators. The first integrator solves for one- dimensional Keplerian orbits using the 
evolution mapping defined by equations flTB]) and the standard Kepler Lagrangian; the sec- 
ond is a general A^-body integrator using the adaptive-stepsize, individual timestep algorithm 
described above. The source code for both integrators is available upon request. 



To illustrate the energy conservation and symplecticity performance of adaptive and 
non-adaptive variational integrators, we performed several one-dimensional simulations of 
the Kepler problem. 



The particular parameters we chose were m = k = 2/9, and / = 2vT^/135, with initial 
conditions tq = 2/45, vq = po = 0.0. With these parameters and initial conditions the total 
energy of the system is = — 1/4, the eccentricity is e = 9/10, and the orbital period is 
T = 2tt (2/45)^''^ ~ 0.0588716. We evolved the system for various times, from r/10 to lOr in 
increments of r/10. The plots shown in Figures H] and [5] use values computed at the end of a 
simulation over the appropriate total time interval, not snapshots of the corresponding values 
in an ongoing simulation; the distinction is important because of what it implies about the 
choices of adaptive timesteps — see Section HI We used three different integration algorithms 
based on equations f|T6|) : a constant timestep integrator, an adaptive timestep integrator 
where h is chosen at the beginning of a step according to 



7. Numerical Experiments 



7.1. 



One-Dimensional Simulation 






(39) 
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where q = r and we set rj = 0.05, and a block-power-of-two timestep integrator which uses 
the above for its /imax- For the constant timestep integrator, we chose the timestep h (ro,Po)- 
All algorithms iterate equation fll6ap to convergence; we thus expect the block-power-of-two 
and constant timestep algorithms to be exactly symplectic. 

Figure H] displays the relative energy error accumulated over many simulated orbits by 
the three algorithms. It is clear that the three algorithms accumulate comparable energy 
error, with the adaptive timestep choice slightly worse than the others. 

Given an evolution mapping Ft : {q,p) ^-> {Qt{q,p), Pt{q,p)), the pushforward of the 
Poincare integral invariant I = dq A dp along the trajectory q{t),p{t) = Ft {qo,Po) is given 
by 

(Ft^) (?o, Po) = [doQt (go, Po) diPt (go, Po) - diQt (go, Po) doPt (go, Po)] dqo A dpo. (40) 



All three integration algor ithms are evolution mappings. An elegant algorithm (jSussman 



20061 : ISussman et al.l l200ll ) exists to compute derivatives of arbitrary computations, such 
as our evolution mappings, without the truncation error which would result from finite 
differencing. We used this algorithm to compute AJ = [{F^J) — J] (go,po) using equation 
(HOj) . A J is a two form, and in a two-dimensional phase space it can be written A J = 
f{t)dqo A dpo, where / is a scalar. The value of / is plotted versus t for the three algorithms 
in Figure O It is clear that the ordinary adaptive stepsize integrator does not conserve 
the symplectic form while the block-power-of-two and constant stepsize integrators do, as 
expected from Section HI As explained in Section |U though the energy error performance 
of the three algorithms is comparable, only the constant-timestep and block-power-of-two 
schemes are symplectic, and only the block-power-of-two scheme is symplectic and allows 
for adaptive timesteps. 



7.2. Many-Body Simulations 



In this subsection we report on the application of the individual and adaptive timestep 
integration algorithm described in this paper to simulations of many-body systems. We 
choose to use the so-called "standard uni ts" : units in which the to tal system mass M = 1, 
G = 1, and the total energy E = —1/4 (jHeggie fc Mathieul 119861 ). In standard units, the 
virial radius of the system is 1. Our initial conditions are a randomly sampled Plummer 
model shifted into a coordinate system where the center of mass is at the origin and the 
total linear momentum is zero. 



Because our code makes no special provision for close encounters, we used a softened 
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Fig. 4. — Relative energy error versus total evolution interval for simulations of the Kepler 
orbit described in the text. The orbital parameters are a = 2/45, e = 9/10, E = —1/4. The 
total evolution interval for each simulation varies from to lOr in steps of r/10. All three 
integration algorithms have comparable energy error. 
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Fig. 5. — Error in the Poincare integral invariant versus total evolution interval for the 
simulations of the Kepler orbit shown in Figure |H The constant timestep (solid line) and 
block-power-of-two timestep (dot-dashed line) algorithms conserve the Poincare integral in- 
variant exactly while the adaptive timestep algorithm (dashed line) doesn't, even though the 
energy error behaviors of the algorithms are comparable. 
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gravitational potential: 



V{t) = ^==, (41) 



with e = 4:/N, where N is the number of bodies in a particular simulation. This is a standard 
techn ique in codes which do not carefully regularize the singular two-body potential (lAarseth 



200l|). 



7.2.1. Symplecticity Tests 

An A^-body system with A^ bodies has a 6A^-dimensional phase space. Given an inte- 
gration mapping, with Jacobian matrix J, conservation of the Poincare integral invariant by 
the mapping implies that 

S = J^SJ, (42) 

where 5* is the "symplectic unit" : 

" -1 
1 



S 



(43) 



where 1 represents the identity matrix in 3A^-dimensions. Using a generalization of the 
algorithm for explicitly computing derivatives of computational algorithms, we can compute 
the Jacobian matrix for the evolution mapping defined by any integration scheme. 

We denote by dl the relative change in the sum of the absolute values of the diagonal 
from lower left to upper right between J^SJ and 5". If equation (H2l) holds, dl = 0. We 
compare dl from the individual and adaptive timestep integrator describe d in this pape r to 



dl from a standard individual and adaptive timestep Hermite integrator (jMakindll99ll ). It 
is difficult to precisely control the stepsizes chosen in an adaptive timestep scheme, so we 
measure dl as a function of the total number of steps taken for the evolution, nsteps- An 
error which scales as h"^ should scale as n^eps- 

Based on the discussion in Section 12.1.2^ we expect that dl from the variational inte- 
grator will be fifth order (that is, it scales as ^^eps); while dl from the Hermite integrator 
will be fourth order (exactly the same order as the integrator's trajectory error). Figure [6] 
demonstrates that this is exactly the case for simulations of a Plummer initial condition with 
N = 25 bodies and varying numbers of steps over a total time interval of T = 1.0 by both 
algorithms. 

Because they both depend on the accuracy of the solution to the implicit equation 
f ll6al) . we expect that symplecticity and angular momentum conservation error would scale 
similarly in a long-term simulation. Results in the next section show that angular momentum 
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dl vs. n . for Hermite and Variational Integrators 
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steps 



Fig. 6. — Error in equation (l42l) (non-conservation of the Poincare integral invariant) in mul- 
tiple integrations of = 25 Plummer initial conditions by the standard Hermite integrator 
and the individual and adaptive timestepping integrator described in the text over a total 
time interval of T = 1.0 with various choices of (adaptive) timestep. The data have been 
least-squares fit to illustrate the scaling with risteps^ the Hermite data fit with dl oc ^^ep^^^, 
and the variational data fit with dl oc ^^ep^s'^^, as expected from the discussion in Section 
[2X21 
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error in a long-running simulation is about one order of magnitude better than the energy 
conservation error (Figure [7]) using our algorithm, while the Hermite algorithm produces 
angular momentum errors commensurate with its energy conservation error. Our algorithm 
appears to have the advantage in angular momentum conservation and symplecticity in 
large- A^, long-time simulations. 



1.2.2. 1000-body Cluster Simulation 
Here we compare a 1000-body cluster simulation run using our varia tional integrator 



with another such simulation using the NB0DY2h code (lAarsethl l200ll ). NB0DY2h is 
not the state-of-the-art in cluster simulations; it is an appropriate comparison to our code 
because it uses the state-of-the-art Hermite integration algorithm but uses softening instead 
of treating close encounters specially. We have run many such simulations — the one described 
here is typical. Recall we use the standard units: the total system mass is M = 1, G = 1, 



and the total energy E = —1/4 (IHeggie fc Mathieul Il986l ). In standard units, the virial 



radius of the system is 1. Initial conditions for both runs are a randomly sampled Plummer 
model shifted into a coordinate system where the center of mass is at the origin and the 
total linear momentum is zero. The random sampling differs between the two codes, so the 
initial conditions are not identical. 

In the simulation reported here, we chose to keep energy conservation to better than a 
part in 5 X 10~^ per unit timestep, with a target of one part in 10^^. If, at the end of an 
advancement by At = 1.0, the relative energy error was greater than 5 x 10~^ we re-started 
the step with a smaller rj parameter (cf. eq. (I57|) ): at the end of every unit time interval, we 
adjusted r] according to 

10-9 

= "aeJe- t^") 

A similar procedure is used to adjust the timestep in NB0DY2h. 

Figure [7] shows the time-evolution of the relative energy error and relative angular 
momentum error in the variational simulation. We can see in Figure [7] that the angular 
momentum conservation error is about an order of magnitude smaller than the energy con- 
servation error. This is consistent with the fact that the angular momentum error scales one 
power of h better than the energy error. 



Figure [H] compares the core radius (computed as described in lCasertano fc HutI (1l985l )) 



versus time for the variational simulation with the core radius output by the NB0DY2h 
simulation (recall that the two simulations do not start with identical initial conditions, but 
rather random samplings of identical phase-space distributions). The dynamical behavior 
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Fig. 7. — Relative energy and angular momentum conservation error for the 1000-body 
simulation performed using the variational integrator. As expected, the angular momentum 
conservation error is much better than the energy conservation error. 
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of the core radius simulated using our algorithm is qualitatively correct. The NB0DY2h 
simulation has momentum (both linear and angular) conservation error which is approx- 
imately commensurate with its energy conservation error — about one order of magnitude 
worse than the variational algorithm's angular momentum conservation error, and six orders 
of magnitude worse than its linear momentum conservation error. Its wall-clock time is sig- 
nificantly better than the varia tional algorithm becaus e our algorithm does not implement 
the nearest-neighbor scheme of I Ahmad fc CohenI (119731 ). and therefore evaluates all 999 po- 
tentials involving the body with the smallest timestep every step, all other 998 potentials 
involving the body with the second- smallest timestep every step for that body, etc. 



8. Discussion 

We have presented a class of A^-body integrators obtained by discretizing the action as 
opposed to discretizing the equations of motion. Such integrators automatically conserve 
discrete momenta and are symplectic. This paper focuses on the fourth-order, three-point 
GGL integrator described in Section I2.1.2[ but Section [2] describes a general framework 
for constructing such integrators. We have demonstrated, theoretically in Section H] and 
experimentally in Section EH that adaptive timestepping integrators can still be symplectic if 
they use a block-power-of-two scheme. Individual timesteps (Section E]) impose requirements 
that reduce the symplecticity and angular momentum conservation from exact to fifth-order 
(one order better than the trajectory error of our integrator); nevertheless, we have shown 
experimentally in Sections 17.2.11 and 17.2.21 that the benefits of fifth-order symplecticity and 
angular momentum conservation are significant. 

Though our code is not CPU-time competitive with standard stellar-cluster simulations 
due to their use of the Ahmad-Cohen nearest-neighbor scheme, we expect it will be useful 
in direct-summation A^-body simulations where phase space volume conservation is more 
important than raw speed — though it would be interesting to see whether the Ahmad- 
Cohen scheme significantly affects the symplecticity of our algorithm in practice, we have 
not investigated this. We intend to use our algorithm to check the accuracy of dark matter 
halo evolution in larger, cosmological A^-body simulations, and to investigate the formation of 
caustics in dark-matter distributions on a smaller scale. In these applications symplecticity is 
essential because the "coldness" of the dark matter phase-space distribution is fundamental 
to the dynamics. 

Additionally, the algorithm was designed with the application to full cosmological A^- 
body simulations in mind. The original motivation for the work was to find a higher- 
order, symplectic algorithm which takes only forward steps as an alternative to the leapfrog 
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Fig. 8. — Core radius versus time for the simulation described in the text using the integrator 
in this paper and an equivalent simulation using the NB0DY2h code. Boxcar averaging has 
been employed with At — 15 to reduce fluctuations in the curves. 



- 30 - 



algorithm used in cosmological simulations. Standard compositional higher-order symplectic 
integrators are not useful in this context because the simulations of gas dynamics in high- 
accuracy cosmological simulations become unstable under evolution backward in time, and 
these integrators all have backwards timesteps. We intend to explore this application in a 
future paper. 

We would hke to thank Jack Wisdom and Gerry Sussman for helpful discussions, encour- 
agement, and inspiration. We would also like to thank an anonymous reviewer for helpful 
comments and for correcting an error in Section [2l This work was supported by NSF grant 
AST-0407050 and NASA grant NNG06-GG99G. 
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